###########################################
# Table A1: Robustness -- ACR with alternative treatment intensity variables

###########################################
options(warn=-1)# removes warnings for perfectly fit observations to remove

dep <- "perc_yes ~"
Xs_elFE <-  paste("turnout","Bond","failed_recently",sep = "+")
FEs <- "|LEAID + election_FE|0"
cluster <- "|LEAID "

dfreg <- df[df$one_election ==0,]
dfreg25 <- dfreg[dfreg$tight_election ==1,]

R0 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,"post:itm_prop_perc_17 +",Xs_elFE,FEs,cluster,sep = "")))
R1 <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,"post:itm_prop_perc_17 +",Xs_elFE,FEs,cluster,sep = "")))
R2 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,"post:ownership_2017 +",Xs_elFE,FEs,cluster,sep = "")))
R3 <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,"post:ownership_2017 +",Xs_elFE,FEs,cluster,sep = "")))

stargazer(R0,R1,R2,R3,
          type = "latex",
          omit.stat = c("ser","f"),
          digits = 2,
          omit.table.layout = "n",
          omit = c("failed_recently","Bond","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA1.alternative_treatments.tex",
          title = "Robustness - ACR with alternative treatment intensity variables",
          label = "tab:alternative_treatment",
          dep.var.labels.include =   FALSE)


rm(list = setdiff(ls(),c("df","city")))

###########################################
# Table A2: Testing for treatment response heterogeneity from SALT cap

###########################################

dfreg <- df[df$one_election==0,]
dfreg25 <- dfreg[dfreg$tight_election==1,]

R1 <- felm(data = dfreg25,weights = dfreg25$population_votingage,
           formula = perc_yes ~ salt_change_ph:post:itm_prop_change + itm_prop_change*post + turnout + Bond + failed_recently | LEAID + election_FE|0|LEAID)
R2 <- felm(data = dfreg25,weights = dfreg25$population_votingage,
           formula = perc_yes ~ salt_amt_change_perc:post:itm_prop_change  + itm_prop_change*post + turnout + Bond + failed_recently | LEAID + election_FE|0|LEAID)
R3 <- felm(data = dfreg25,weights = dfreg25$population_votingage,
           formula = perc_yes ~ wasted_salt_ph:post:itm_prop_change  + itm_prop_change*post + turnout + Bond + failed_recently | LEAID + election_FE|0|LEAID)
R4 <- felm(data = dfreg25,weights = dfreg25$population_votingage,
           formula = perc_yes ~ wasted_salt_perc:post:itm_prop_change  + itm_prop_change*post + turnout + Bond + failed_recently | LEAID + election_FE|0|LEAID)



stargazer(R1,R2,R3,R4,type = "latex",
          omit.stat = c("ser","f"),
          omit = c("failed_recently","Bond","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA2.robust_SALT.tex",
          title = "Testing for treatment response heterogeneity from SALT cap",
          label = "tab:refdum_results_intensive",
          digits = 2)


rm(list = setdiff(ls(),c("df","city")))

###########################################
# Table A3: Placebo regressions testing the “Yes” votes share parallel trend

###########################################

# Placebo by changing the year of the treatment
placebo <- df[df$year<=2018,]
a <- as.data.frame(table(placebo$LEAID))
a <- a[a$Freq>1,]$Var1
placebo <- placebo[placebo$LEAID %in% a,]

count <- 1
R <- list()
for (i in c(2009,2010,2011,2012,2013,2014,2015,2016,2017,2018)){
  
  placebo$post <- ifelse(placebo$year >= i,1,0)
  
  placebo_reg <- placebo[placebo$tight_election ==1,]
  
  R[[count]] <- felm(data = placebo_reg,weights = placebo_reg$population_votingage,
                     formula = perc_yes ~ itm_prop_change*post + turnout + Bond + failed_recently | LEAID + election_FE|0|LEAID)
  
  count <- count +1
}


stargazer(R,type = "latex",
          omit.stat = c("ser","f"),
          column.labels = as.character(seq(2009,2018,1)),
          omit = c("failed_recently","Bond","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA3.placebo_year.tex",
          title = "Placebo regressions testing the ``Yes'' votes share parallel trend",
          label = "tab:placebo_year",
          digits = 2)



rm(list = setdiff(ls(),c("df","city")))

###########################################
# Table A4: Robustness - Covid-19 and virtual mode of instruction

###########################################

## Merge with mode of instructions in schools

sch.mode <- read.csv("Data/Raw data/Instructional_modes.csv")

sch.mode <- sch.mode[sch.mode$StateAbbrev == "CA",!(names(sch.mode)%in% c("DistrictName","StateAssignedDistrictID","StateAbbrev"))]
colnames(sch.mode)[1] <- "LEAID"

## Some school districts have different LEAID in the mode of instruction database
sch.mode$LEAID[which(sch.mode$LEAID == 601426)] <- 605670# Bonsall school district
sch.mode$LEAID[which(sch.mode$LEAID == 600153)] <- 606037# Alhambra school district
sch.mode$LEAID[which(sch.mode$LEAID == 601442)] <- 640500# Upper Lake school district
sch.mode$LEAID[which(sch.mode$LEAID == 601428)] <- 642990# Wiseburn school district


sch.mode$LEAID <- sprintf("%07d",sch.mode$LEAID)

df2 <- merge(df,sch.mode,by = "LEAID",all.x = TRUE)

## Regressions

dep <- "perc_yes ~"
Xs_elFE <-  paste("turnout","Bond","failed_recently",sep = "+")
FEs <- "|LEAID + election_FE|0"
cluster <- "|LEAID "


dfreg <- df2[df2$one_election==0,]
dfreg25 <- dfreg[dfreg$tight_election==1,]

main <- "post*share_virtual +"
R0 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R1 <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))


main <- "post*share_virtual + post*itm_prop_change +"
R2 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R3 <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))


stargazer(R0,R1,R2,R3,
          type = "latex",
          omit.stat = c("ser","f"),
          digits = 2,
          omit.table.layout = "n",
          report = "vc*s",
          omit = c("failed_recently","Bond","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA4.robust_school_closure.tex",
          title = "Robustness - Covid-19 and virtual mode of instruction",
          label = "tab:virtual_instructions",
          dep.var.labels.include =   FALSE)

rm(list = setdiff(ls(),c("df","city")))

###########################################
# Table A5: Voters’ tax deductions and local governments’ margin of adjustment

###########################################

## PART I OF THE TABLE A5

source("Rscripts/11bis.construct_matrix_of_refdums.R")

R_pois <- glm(nb_refdum ~ as.factor(LEAID) + as.factor(year) + itm_prop_change:post, data=rdm, family ="poisson")
R_bin <- glm(refdum_binary~as.factor(LEAID) + as.factor(year) + itm_prop_change:post,data=rdm, binomial(link = "probit"))


## PART II OF THE TABLE A5
dfreg <- df[df$one_election==0,]

dfreg$bond_ph <- dfreg$bond_ph/1000 #For readibility
dfreg$levy_parcel <- dfreg$levy_parcel/100#For readibility


R0 <- felm(data = dfreg,weights = dfreg$population_votingage,
           formula = bond_ph ~ itm_prop_change:post  | LEAID + election_FE|0|LEAID)
R1 <- felm(data = dfreg,weights = dfreg$population_votingage,
              formula = levy_increase ~ itm_prop_change:post  | LEAID + election_FE|0|LEAID)
R2 <- felm(data = dfreg,weights = dfreg$population_votingage,
           formula = turnout ~ itm_prop_change:post  | LEAID + election_FE|0|LEAID)


stargazer(R_bin,R_pois,R2,R0,R1,
          type = "latex",
          omit.stat = c("ser","f","aic"),
          omit = c("Constant","year","LEAID"),
          out = "Output/TableA5.balance_test.tex",
          title = "Voters' tax deductions and local governments' margin of adjustment",
          label = "tab:balance_tests",
          digits = 2)

## To get marginal effects of the two first model (as shown in the appendix), run following model
  #library(mfx) # to compute marginal effect in binary/poisson models
  #a <- logitmfx(formula = refdum_binary~as.factor(LEAID) + as.factor(year) + itm_prop_change:post,data =rdm)
  #a$mfxest['itm_prop_change:post',]
  #b <- poissonmfx(formula = nb_refdum ~ as.factor(LEAID) + as.factor(year) + itm_prop_change:post,data = rdm)
  #b$mfxest['itm_prop_change:post',]


rm(list = setdiff(ls(),c("df","city")))


###########################################
# Table A6: Average Causal Response by Frequency of Referendums

###########################################

dep <- "perc_yes ~"
main <- "post:itm_prop_change +"
Xs_elFE <-  paste("turnout","Bond","failed_recently",sep = "+")
FEs <- "|LEAID + election_FE|0"
cluster <- "|LEAID "

a <- as.data.frame(table(df$LEAID))

dfreg0 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=10,]$Var1,]
dfreg1 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=9,]$Var1,]
dfreg2 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=8,]$Var1,]
dfreg3 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=7,]$Var1,]
dfreg4 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=6,]$Var1,]
dfreg5 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=5,]$Var1,]
dfreg6 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=4,]$Var1,]
dfreg7 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=3,]$Var1,]
dfreg8 <- df[df$LEAID %in% a[a$Freq>1 & a$Freq <=2,]$Var1,]

R0 <- felm(data = dfreg0,weights = dfreg0$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R1 <- felm(data = dfreg1,weights = dfreg1$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R2 <- felm(data = dfreg2,weights = dfreg2$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R3 <- felm(data = dfreg3,weights = dfreg3$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R4 <- felm(data = dfreg4,weights = dfreg4$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R5 <- felm(data = dfreg5,weights = dfreg5$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R6 <- felm(data = dfreg6,weights = dfreg6$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R7 <- felm(data = dfreg7,weights = dfreg7$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R8 <- felm(data = dfreg8,weights = dfreg8$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))


stargazer(R0,R1,R2,R3,R4,R5,R6,R7,R8,
          type = "latex",
          omit.stat = c("ser","f"),
          digits = 2,
          omit.table.layout = "n",
          omit = c("failed_recently","Bond","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA6.results_by_frequency.tex",
          title = "Average Causal Response by Frequency of Referendums",
          label = "tab:placebo_frequency",
          dep.var.labels.include =   FALSE)


rm(list = setdiff(ls(),c("df","city")))


###########################################
# Table A7: City Referendums - Placebo test

###########################################


dep <- "percent_yes ~"
main <- "chg_prop.itm*post +"
Xs <- paste("presidential_el","midterm_el","odd_year","turnout","failed_recently","type",sep = "+")
Xs_elFE <-  paste("turnout","type","failed_recently",sep = "+")
Xs_sales <-  paste("turnout","failed_recently","SalesTax",sep = "+")
FEs <- "|City_geoid + election_FE|0"
cluster <- "|City_geoid  "

# Removes unique cities and create close elections samples
a <- as.data.frame(table(city$City_geoid))
a <- a[a$Freq>1,]
dfreg <- city[city$City_geoid %in% a$Var1,]
rm(a)
dfreg25 <- dfreg[dfreg$percent_yes >=40 & dfreg$percent_yes<=90,]

R0 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,"post+" ,Xs,"|City_geoid |0",cluster,sep = "")))
R1 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,main,Xs,"|City_geoid |0",sep = "")))
R2 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R3 <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))
R4 <- felm(data = dfreg25[dfreg25$type == "Sales Tax" ,],weights = dfreg25[dfreg25$type == "Sales Tax" ,]$population_votingage,formula = as.formula(paste(dep,main,Xs_sales,FEs,cluster,sep = "")))

stargazer(R0,R1,R2,R3,R4,
          type = "latex",
          omit.stat = c("ser","f"),
          digits = 2,
          omit.table.layout = "n",
          omit = c("failed_recently","type","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableA7.placebo_city.tex",
          title = "City Referendums - Placebo test",
          label = "tab:refdum_city",
          dep.var.labels.include =   FALSE)


rm(list = setdiff(ls(),c("df","city")))


###########################################
# Table B1: Robustness results for bond only referendums

###########################################


dfreg <- df[df$one_election ==0,]
dfreg25 <- dfreg[dfreg$tight_election ==1,]

dep <- "perc_yes ~"
main <- "post:itm_prop_change +"
FEs <- "|LEAID + election_FE|0"
cluster <- "|LEAID "

# Excludes control for bond amount
Xs_bond <-  paste("turnout","failed_recently",sep = "+")
R0 <- felm(data = dfreg[dfreg$type == "Bond" ,],weights = dfreg[dfreg$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep,main,Xs_bond,FEs,cluster,sep = "")))
R1 <- felm(data = dfreg25[dfreg25$type == "Bond" ,],weights = dfreg25[dfreg25$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep,main,Xs_bond,FEs,cluster,sep = "")))

# Includes control for bond amount in log form
Xs_bond <-  paste("turnout","log(bond_ph)","failed_recently",sep = "+")
R2 <- felm(data = dfreg[dfreg$type == "Bond" ,],weights = dfreg[dfreg$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep, main,Xs_bond,FEs,cluster,sep = "")))
R3 <- felm(data = dfreg25[dfreg25$type == "Bond" ,],weights = dfreg25[dfreg25$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep, main,Xs_bond,FEs,cluster,sep = "")))

# Includes control for bond amount 
Xs_bond <-  paste("turnout","bond_ph","failed_recently",sep = "+")
R4 <- felm(data = dfreg[dfreg$type == "Bond" ,],weights = dfreg[dfreg$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep,main,Xs_bond,FEs,cluster,sep = "")))
R5 <- felm(data = dfreg25[dfreg25$type == "Bond" ,],weights = dfreg25[dfreg25$type == "Bond" ,]$population_votingage,formula = as.formula(paste(dep,main,Xs_bond,FEs,cluster,sep = "")))

stargazer(R0,R1,R2,R3,R4,R5,
          type = "latex",
          omit.stat = c("ser","f"),
          digits = 2,
          omit.table.layout = "n",
          omit = c("failed_recently","type","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableB1.robust_bondonly.tex",
          title = "Robustness results for bond only referendums",
          label = "tab:bond_only",
          dep.var.labels.include =   FALSE)


rm(list = setdiff(ls(),c("df","city")))

###########################################
# Table B2: Differences in approval rates for bonds and other referendums – Testing the effects of housing supply elasticity

###########################################


land <- read.csv("Data/Raw data/land_use_school_districts_all.csv")
land$NCESID <- sprintf("%07d",land$NCESID)
land$state <- substr(land$NCESID,1,2)
land <- land[land$state =="06",]
land$developable_perc <- 1-land$developed_perc
land <- land[,c("NCESID","developed_perc","developable_perc")]

df_with_land <- merge(df,land,by.x = "LEAID",by.y = "NCESID")
rm(land)

dfreg <- df_with_land[df_with_land$one_election ==0 & df_with_land$post ==0,]

R0 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = perc_yes ~ Bond*developable_perc + presidential_el + midterm_el + odd_year + turnout + failed_recently|0|0|LEAID )

R1 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = perc_yes ~ Bond*developable_perc + turnout + failed_recently| election_FE |0|LEAID)

R2 <- felm(data = dfreg,weights = dfreg$population_votingage,formula = perc_yes ~ Bond*developable_perc + turnout + failed_recently| election_FE + LEAID|0|LEAID)

stargazer(R0,R1,R2,
          type = "latex",
          omit.stat = c("f","ser"),
          omit = c("failed_recently","type","presidential_el","midterm_el","odd_year","turnout"),
          out = "Output/TableB2.bond_with_housingsupply.tex",
          title = "Differences in approval rates for bonds and other referendums -- Testing the effects of housing supply elasticity",
          label = "tab:bond_housing_supply",
          digits = 2)


rm(list = setdiff(ls(),c("df","city")))

options(warn=0)

